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The problem of computing a representation for a real polynomial as a sum of minimum number 
of squares of polynomials can be casted as finding a symmetric positive semidefinite real matrix 
(Gram matrix) of minimum rank subject to linear equality constraints. In this paper, we propose 
algorithms for solving the minimum-rank Gram matrix completion problem, and show the con- 
vergence of these algorithms. Our methods are based on the modified fixed point continuation 
(FPC) method. We also use the Barzilai-Borwein (BB) technique and a specific linear combi- 
nation of two previous iterates to accelerate the convergence of modified FPC algorithms. We 
demonstrate the effectiveness of our algorithms for computing approximate and exact rational 
sum of squares (SOS) decompositions of polynomials with rational coefficients. 

1. Introduction 

Let x = [xi, . . . , x s ] and / 6 RLr], then / is a sum of squares (SOS) in R[x] if and only if it can be 
written in the form 



in which rrid(x) is a column vector of monomials of degree less than or equal to d and W is a real 
positive semidefinite matrix [39 \ Theorem 1] (see also pZ3]). W is also called a Gram matrix for /. 
If W has rational entries, then / is a sum of squares in Qfaq, . . . ,x n ]. 

Problem 1. Let f G Q[xi, • • • ,x s ] be a polynomial of the degree 2d, compute a representation for 
it as a sum of minimum number of squares of polynomials in Q[x\, . . . ,x s ]. 

The set of all matrices W for which (JT]) holds is an afflne subspace of the set of symmetric 
matrices. If the intersection of this affine subspace with the cone of positive semidefinite (PSD) 
matrices is nonempty, then / can be written as a sum of squares. Since the components of m^x) 
are not algebraically independent, W is in general not unique. Problem Q] can be restated as finding 
a Gram matrix with minimum rank satisfying a given set of constraints: 



For s = 1, Pourchet's main theorem [38] implies that every positive definite univariate polyno- 
mial in Q[x] is a sum of five squares in Q[x]. Therefore, the minimum rank of the Gram matrix 
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f(x) = m d (x) ■ W ■ m d (x) 



(1) 
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satisfying ([2]) is bounded by 5 for s = 1. For s > 1, Pfister's general theorem [37] shows that 
every positive definite polynomial in M.[x±, . . . ,x s ] is a sum of 2 s squares of rational functions in 
. . . ,x s ). It is well known that there exist positive semidefinite polynomials which cannot be 
written as sums of polynomial squares. However, as shown in [21 j . various exceptional SOS prob- 
lems in the literature by Motzkin, Delzell, Reznick, Leep and Starr, the IMO'71 problem by A. 
Lax and P. Lax, and the polynomial Vor2 in can be written as sums of less than 10 squares 
of polynomials after multiplying by suitable polynomials. The advantage of computing a numerical 
Gram matrix with small rank is that we can refine the approximately computed Gram matrix to 
high accuracy by structure preserved Gauss-Newton iteration more efficiently in order to recover 
the exact SOS representation of / [20" 1 |2~T [ l35 | [36] . 

In general, the rank minimization is an intractable problem and is in fact provably NP-hard due 
to the combinational nature of the non-convex rank function [9] . In [121 [TBI [41] , they showed that 
rank(W) can be replaced by the nuclear norm of W, which is the best convex approximation of the 
rank function over the unit ball of matrices. Expanding the right-hand side of the equality condition 
of ([21), matching coefficients of the monomials, we obtain a set of linear equations for the entries 
of W which can be written as .A(W) = b, where the action of the linear operator A : S n — > W on 
W is described by Tv(AjW),i = 1, . . . ,p for A u . . . , A p e R nxn . We use A* : W ->• S n to denote 
the adjoint operator of A. The rank minimization problem ([2]) can be relaxed to the nuclear norm 
minimization problem 



min ||W||* 
s.t. A(W) = b 

w y o, w T = w 



(3) 



where the nuclear norm || W||* is defined as the sum of its singular values. The constraint «4(W) = b 
can also be relaxed, resulting in either the problem 

min ||W||* 

s.t. \\A(W) - b\\ 2 < e } (4) 



or its Lagrangian version 



w y 0, W T = W 



min ^\\W\U + \\\A{W)-b\\l (5) 



where S" is the set of symmetric positive semidefinite matrices and > is a parameter. 

In [H [161 [Ml [23 [26] , they studied how to determine whether partially specified positive semidef- 
inite matrices can be completed to fully specified matrices satisfying certain prescribed properties. 
A number of recent work has also shown that the low-rank solution can be recovered exactly via 
minimizing the nuclear norm under certain conditions [7J El BTl W2\ . Several algorithms based on 
the interior point method have been proposed in [H [28| HT| |4*3"1 |4"5] for solving the semidefi- 
nite programming problem derived from the rank minimization problem Q. Since most of these 
methods use second-order information, the memory requirement for computing descent directions 
quickly becomes too large as the problem size increases. Recently, several fast algorithms using only 
first-order information have been developed in [61 [T4"l [T9| [30] [3T| 147] . These first-order methods, 
based on function values and gradient evaluation, cannot yield as high accuracy as interior point 
methods, but much larger problems can be solved since no second-order information needs to be 
computed and stored. 

Motivated by these exciting work, in this paper, we present two algorithms for solving the 
minimum-rank Gram matrix completion problem ([3]). Our algorithms are based on the modified 
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fixed point continuation method. By modifying the shrinkage operator in FPC and using the 
Barzilai-Borwein technique to compute explicit dynamically updated step sizes, we get an algorithm, 
called modified fixed point continuation method with the Barzilai-Borwein technique (MFPC-BB). 
We prove the convergence of our algorithm under certain condition. Some accelerated gradient 
algorithms were proposed in [U [191 [32l [33j [Ml [47J [49] . These algorithms rely on computing the next 
iterate based not only on the previous one, but also on two or more previously computed iterates. 
These accelerated gradient methods have an attractive convergence rate of 0(1/ k 2 ), where k is the 
iteration counter. We incorporate this accelerating technique in the MFPC-BB algorithm to get 
an accelerated fixed point continuation algorithm with the Barzilai-Borwein technique (AFPC-BB), 
which shares the improved rate 0(1/ k 2 ) of the optimal gradient method. 

We also notice that algorithms in the literature mostly focus on recovering a randomly generated 
large-scale matrix from incomplete samples of its entries. Although it has been pointed out briefly 
in |47j that these algorithms can be adapted easily to solve the regularized semidefinite linear 
least squares problem ([5]), it is interesting for us to investigate how to use these newly developed 
techniques to compute approximate and exact rational sum of squares (SOS) decompositions of 
polynomials with rational coefficients. 

Notations: Let § n C M nxn denote the space of symmetric n x n matrices. The inner product 
between two elements X, Y 6 § n is denoted by (X, Y) = Tv(X T Y). The Frobenius norm of a matrix 
X is denoted by ||^||f, the nuclear norm by \\X\\* and the operator norm (or spectral norm) by 

The rest of the paper is organized as follows. In Section 2, we derive the modified fixed point 
iterative algorithm for the minimum-rank Gram matrix completion problem. In Section 3, we 
establish the convergence result for the iterations given in Section 2 and prove that it converges to 
the optimal solution of the regularized linear least squares problem ([5]). In Section 4, we introduce 
two techniques to accelerate the convergence of our algorithm and present MFPC-BB and AFPC- 
BB algorithms for solving problem ([5]). We demonstrate the performance and effectiveness of 
our algorithms through numerical examples for computing approximate and exact rational sum of 
squares decompositions of polynomials with rational coefficients in Section 5. 

2. Modified fixed point iterative algorithm 

Let / : M n i xn 2 — > R be a convex function, the subdifferential of / at X* G R"i xri 2 denoted by df 
is the compact convex set defined by 

df(X*) := {Z E R™i x ™2 : f(Y) > f{X*) + (Z,Y-X*),V Y € M niX " 2 }. 

Following discussions in \27\ Theorem 3.1] and [51j . we derive the expression of the subdifferential 
of the nuclear norm at a symmetric matrix. 

Theorem 1. Let W G then 

d\\W\\* = {Q [1) Q {1)T - Q (2) Q (2)T + Z : Q^ T Z = 0,i = 1, 2, and \\Z\\ 2 < 1}, 

where and are orthogonal eigenvectors associated with the positive and negative eigenvalues 
of W respectively. 

Proof. Suppose that the eigenvalues of a symmetric matrix W can be ordered as Ai > • • • > At > 
> At + i > • • • > A s , A s+ i = • • • = A n = 0. Let W = QAQ T be a Schur decomposition of W, where 
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Q £ M nxn is an orthogonal matrix and A = diag(Ai, . . . , A n ). These matrices can be partitioned as 

/ A« \ 

Q=(Q {1) ,Q {2 \Q {3) ), a= o A( 2 ) , (6) 

V ' \ A( 3 ) / 

with Q^ 2 \ having t, s — t,n — s columns and being associated with A^ = diag(Ai, . . . , X t ), 
A^ 2 ) = diag(At+i, . . . , A s ), and A^ 3 ) = diag(A s +i, . . . , A„), respectively. 
Let A = (Ai, . . . , X n ) T and recall that 

d\\\\\i = {y G M n : yi = l,i = l,...,t; Vj = =t + l,...,s; \y k \ < 1, k = s + 1, . . . , n}. 
Let Y" G by [271 Theorem 3.1], we have 

Y = Qdiag(d) Q T , 

where d S 9||A[|i. Therefore 

Y = QWQW - Q<?)Q&) T + Q^DQW, 

where D is an (n — s) x (n — s) diagonal matrix with diagonal elements less than 1 in modulus. 

Let Z = QWDQWr, we have Q^ T Z = 0,i = 1,2. Let cxi(-) denote the largest singular value 
of a given matrix, then we have 

\\Zh = Q {3) DQ^ T < ai (D) < 1, 

which completes the proof. □ 
The optimality condition in [3Q|, Theorem 2] can be generalized to the optimality condition for 
the constrained convex optimization problem ([5]). 

Theorem 2. Let f : S n — > R be a proper convex function, i.e. f < +oo for at least one point and 
f > — oo for every point in its domain. Then W* is an optimal solution to the problem 

min f(W) (7) 

if and only ifW* £ S™, and there exists a matrix U £ df(W*) such that 

(U, V - W*) > 0, for all V € §+. (8) 
Proof. Suppose U € df(W*) and satisfies the inequality condition (jHJ), hence 

f(v) > f(w*) + (u,v- w*), vye§;, 

we have f(V) > f(W*), for all V £ S™ . This shows that VF* is an optimal solution of the problem 
©• 

Conversely, suppose W* is the optimal solution of the problem (J7J), and flSJ) does not hold, i.e., 
there exists f7 G 9/(W*), such that 

3 VgS", s.t. (U, V - W*) < 0. (9) 

Consider = tW* + (1 — t)V, where t G [0, 1] is a parameter. Since Z(t) is on the line segment 
between W* and V, and S2. is a convex set, Z(i) G S",V i G [0,1]. By @U Theorem 23.4], the 
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one-sided directional derivative of / at Z(l) with respect to the vector W* — V satisfies the following 
equation 

f(Z(t); W* - V)\ t=1 = f'(W*; W*-V)= sup{{W, W* — V) : W £ df{W*)}. 
According to ([9]), we have 

f'(Z(t); W* - V)\ t=l > {U, W* — V}> 0. 

Therefore, for a small value e > 0, we have f(Z(l — e)) < f(W*), which is contradict to the fact 
that W* is optimal to the problem ([7]). □ 
Based on above theorems, we can introduce a thresholding operator and extend the fixed point 
iterative scheme for solving ([5]). 

Definition 1. Suppose W = QAQ T is a Schur decomposition of a matrix W S § n , where A = 
diag(Ai, . . . , A n ) and Q is a real orthogonal matrix. For any v > 0, the matrix thresholding operator 
%,{■) is defined as 

%(W) := Q%(A)Q T , %(A) = diag({Ai - v}+), 

where t+ = max(0,t). 

We should point out that the idea of using the eigenvalue decomposition of Y k has also appeared 
in [47} Remark 3]. However, to our best knowledge, there exists no convergence analysis about the 
eigenvalue thresholding operator in the literature. 

Let fi and r be positive real numbers and X° be an initial starting matrix. For k = 0, 1, 2, • • • , 
we compute 



Y k = X k -rA*(A(X k ) -b), 
until a stopping criterion is reached. 



X k+l = T~ Tfl (Y k ), ' ' ' (10) 



Theorem 3. Suppose a matrix W* £ S™ satisfies 

1. \\AiyV*) — &H2 < lijn for a small positive number \x. 

2. W* = T T ^{h{W*)), where h(-) = /(•) - tA*(A(-) - b) and /(•) is an identity operator. 

Then W* is the unique optimal solution of the problem |5|j. 

Proof. Let v = r\i and Y* = h(W*) = W* + E G S n , where E = -tA*(A(W*) - b). We claim 
that T V (Y*) i s the unique optimal solution to the following problem 

min i/||W||* + -||W-y*|||. > (11) 
w&i 2 

In fact, since the objective function ^||W||* + ^\\W — Y*\\ 2 F is strictly convex, there exists a 
unique minimizer, and we only need to prove that it is equal to T V (Y*). Without loss of generality, 
we assume that the eigenvalues of Y* can be ordered as 

AiOn > > X t (Y*) >u> A t+1 (Y*) >...>()>•••> X S (Y*), X S+1 (Y*) = ■■■ = X n (Y*) = 0. 
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We compute a Schur decomposition of Y* as 

y* = qWaWqW t + qWaWqW t , 

where A« = diag(Ai, . . . , At), A (2) = diag(A m , . . . , A s ), Q« and are block matrices corre- 
sponding to A^ 1 ) and A^ 2 ) respectively. Let X = T U {Y*), we have 

X = Q^\A^-uI)Q^ T , 

therefore, 

Y* - X = v(qWqW t + Z), Z = v- 1 Q&>A®QW r . 

By definition, QW t Z = 0. 

• If \ t+1 (Y*) > \X S (Y*)\, then \\Z\\ 2 = X t+ i(Y*)/u < 1. 

• Otherwise, let y = (y 1} . . . , y p ) T = A(W*) -beW, then 

II 7tm|2 2 1 1 a* i|2 ^ 2 2/ 2 , , 2\ ^ 2 2 

\E\ F = t \\A y\\ F < t n (y 1 H h y p ) < r // . 

Notice that £ G S™ and G S£, by [15, Theorem 8.1.5], we have 

= |A 8 (y*)| = maxjlA^^UAn^)!} < JjSjj^ < ± 

V V ~ V 

Hence, according to Theorem[H we have Y* —X G which means that G — y*. 

By Theorem [21 we immediately conclude that T V {Y*) is an optimal solution of the problem (jlip . 

Since the objective function of the problem ([5]) is strictly convex, its optimal solution is also 
unique. If W* = %n(Y*), by Theorem [21 there exists a matrix U G + W* - Y* such that 

(u, v - w*) > o, v y g s+. 

Let U = U/t, by substituting z/ = r/i and y* = W* — 7v4* (.A(W*) -6) into the above subdifferential 
function, we have U G //0||W*[|* + A*(A(W*) - b) satisfying 

{u, v - w*) > o, vyG§+. 

By applying Theorem [2] once again, it is true that W* is the optimal solution of the problem ©. □ 

3. Convergence analysis 

In this section, we analyze the convergence properties of the modified fixed point iterative scheme 
(llOp . We begin by recording two lemmas which establish the non-expansivity of the thresholding 
operator 7~u(h(-)). 

Lemma 1. The thresholding operator T v is non-expansive, i.e., for any X\,X 2 G W 1 , 

\\%{Xi) - %(X 2 )\\ F < - X 2 \\ F . (12) 

Moreover, 

||*i - X 2 \\ F = \\%(X 1 ) - %(X 2 )\\ F ^X t -X 2 = %(Xx) - %(X 2 ). 
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Proof. Let X 1 = QWaWqW t and X 2 = Q^A^Q^ T be Schur decompositions of X 1 and X 2 , 
respectively, where 

{1) _ / diag(Ai) \ ( 2 ) _ / diag(A 2 ) 
V° / V 

Ai = (ai, . . . ,a s ) T and A2 = (f3\, . . . , (3t) T are vectors of eigenvalues of X\ and X2 respectively, 
and , are orthogonal matrices. Suppose that a.\ > ■ ■ ■ > > v > afc+i > • • • > a s and 
/?i > • • • > A > v > > ■ ■ ■ > Pti then we have 

X X := %(Xi) = Q^A^QW, X 2 := %{X 2 ) = Q^A^Q^ 7 \ 

where 

7(1) = ( diag(Ai) A ^ (2) = ( diag(A 2 ) 
V J ' V 

Ai = («i — v, . . . , otk — v) and A 2 = — v,. . . ,f3\ — u). Therefore, we have 

\\x l -x 2 f F -\\x l -x 2 \\ 2 F 

= Tr((Xx - X 2 f(X l - X 2 )) - Tr((Xi - X 2 f \X t - X 2 )) 

= Tr(XfX! - XfX! + XjX 2 - XjX 2 ) - 2Ti(XfX 2 - XfX 2 ) 

= B<* " -) 2 + E ^ " ECA - -) 2 - 2Tr(XfX 2 - XfX 2 ). 

i=i i=i i=i i=i 

It is known that for symmetric matrices X, Y, 

Tv(XY) < \(X) T X(Y), 

with equality if and only if there exists an orthogonal matrix Q such that 

X = Qdiag(A(X))Q T , Y = Qdiag(A(y))Q T , 

where X(X),X(Y) are the vectors of eigenvalues of X and Y respectively (see j27J Theorem 2.2]). 
Hence, without loss of generality, assuming k < I < s < t, we have 

Tr(XfX 2 - XfX 2 ) = Tr((Xx - X 1 ) T (X 2 - X 2 ) + (X 1 - X X ) T X 2 + Xf(X 2 - X 2 )) 

< \(X! - Xi) T X{X 2 - X 2 ) + X{X 1 - X!) T A(X 2 ) + X(X 1 ) T X(X 2 - X 2 ) 

I s k I 

1=1 i=l+l i=l i=k+l 

Therefore, 



\\x x - x 2 \\i - \\x x - x 2 \\i > «? - £(«* - v ) 2 + E P - £(& 

i=l 8=1 i=l i=l 

I s k I 

-2E«i^+ E «ift + £(ft-^ + E 
i=i «=;+i i=l i=fc+i 

s t s £ 

^ ( E « 2 + E $ - 2 E ai &) + E ^ - v2 + ° 2 - 2a ^*)- 



i=Z+l i=i+l i=2+l i=fc+l 
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Since t > s and a 2 + j3 2 — 2aifa > 0, we obtain 

s t s 

E a ? + E & - 2 E ^ °- 

i=2+l i=Z+l i=2+l 

Moreover, since the function g(x) = 2fax — x 2 + cq — 2otifa is monotonically increasing in [— oo, fa], 
and a.i < v < fa, i = k + 1, . . . , I, 

2vfa - v 2 + a 2 - 2aifa > 0, i = k + l,...,l. 

Hence, we have 

\\Xl ~~ ^2||f ~~ ||-X~1 ~~ ^2||f — 0' 

i.e., (H2D holds. 

Furthermore, if ||-Xi — X 2 \\ F = \\T u (Xi) — T u (Xi)\\f, then s = t, k = I and ai = fa,i = fe+1, . . . , s, 
which further implies that AW - A.M = A^ 2 ) - A^ 2 ) and Tr((Xi - X 1 ) T (X 2 - X 2 )) achieves its 
maximum. Hence, there exists an orthogonal matrix Q such that 

X 1 -X 1 = Q(A« - A«)Q T = Q(A^ - 1^)Q T =X 2 - X 2 , 

which implies that 

X 1 -X 2 = %(X 1 )-%(X 2 ). (13) 

Suppose (fT3|) holds, then \\Xi — X 2 \\ F = \\T u (Xi) — T U {X 2 ) which completes the proof. □ 
The following lemma and its proof are analogous to results in \17\ [30] . 

Lemma 2. Suppose that the step size r satisfies r £ (0, 2/||„4.|| 2 ). Then the operator h(-) = /(•) — 
tA*(A(-) — b) is non-expansive, i.e., for any X\,X 2 £ § n ; 

\\h{X x )-h{X 2 )\\ F <\\X x -X 2 \\ F . 

Moreover, we have 

\\h(X!) - h(X 2 )\\ F = \\X X - X 2 \\ F h(X{) - h(X 2 ) = x 1 - x 2 , 

where /(•) is an identity operator. 

We now claim that the modified fixed point iterations ()10p converge to the optimal solution of 
the problem ([5]). 

Theorem 4. Let r £ (0, 2/||^4|||) and W* £ satisfy 

1. \\A{W*) — b\\ 2 < ji/n for a small positive number \i. 

2. W* = r Tli (h(W*)), where h(-) = !(•) - tA*(A(-) - b). 

Then the sequence {X k } obtained via modified fixed point iterations J7Q|) converges to W*. 
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Proof. Let v = r\i. Since both T v {-) and h(-) are non-expansive, T u (h(-)) is also non-expansive. 

Therefore, {A fc } lies in a compact set and must have a limit point. Suppose X = linij >00 X fc J 

satisfying \\A(X) - b\\ 2 < /u/n. By W* = %(h{W*)), we have 

\\X k+1 - W*\\ F = \\%(h(X k )) - %{h{W*))\\ F < \\h(X k ) - h(W*)\\ F < \\X k - W*\\ F , 

which means that the sequence {||A fc — W*||.f} is monotonically non-increasing. Therefore 

lim \\X k - W*\\ F = \\X - W*\\ F , 

k — >oo 

where X can be any limit point of {X k }. By the continuity of 7^(/i(-))> we have 

%(h{X))= lim %(h(X k >))= lim X k ? + \ 

j — >oo j — >oo 

i.e., T v (h{X)) is also a limit point of {X k }. Therefore, we have 

\\%(h(X)) - %(h(W*))\\ F = \\%(h(X)) - W*\\ F = \\x - w*\\ F . 
Using Lemma Q] and Lemma [2] we obtain 

%(h(X)) - %(h(W*)) = h(X) - h{W*) =X-W*, 

which implies T u (h(X)) = X. By Theorem El X is the optimal solution to the problem ([5]), i.e., 
X = W*. Hence, we have 

lim \\X k -W*\\f = 0, 

k — >-oo 

i.e., converges to its unique limit point W*. □ 

4. Implementation 

This section provides implementation details of the modified FPC algorithm for solving the minimum- 
rank Gram matrix completion problem. 

4.1. Evaluation of the eigenvalue thresholding operator 

The main computational cost of the modified FPC algorithm is computing the Schur decompositions. 
Following the strategies in [6JH7], we use PRO PACK [23] in Matlab to compute a partial Schur 
decomposition of a symmetric matrix. 

PROPACK can not automatically compute only eigenvalues greater than a given threshold v. 
To use this package, we must predetermine the number Sk of eigenvalues of Y k to compute at the 
/c-th iteration. Suppose X k = Q k ~ l A fc_1 (Q k ~ 1 ) T , we set equal to the number of diagonal entries 
of A fc_1 that are no less than £fc||A fc-1 ||2, where is a small positive number. Notice that s/% is 
non-increasing. If s& is too small, the non-expansive property (1120 of the thresholding operator T v 
may be violated. We increase s/% by 1 if the non-expansive property is violated 10 times |30j. 
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4.2. Barzilai-Borwein technique 

In [30], the authors always set the parameter r = 1 since their operator A is generated by randomly 
sampling a subset of p entries from matrices with i.i.d. standard Gaussian entries. For this linear 
map, the Lipschitz constant for the objective function of © is 1. According to Theorem HI conver- 
gence for the Gram matrix completion problem is guaranteed provided that r 6 (0, 2/||^4||2). This 
choice is, however, too conservative and the convergence is typically slow. 

There are many ways to select a step size. For simplicity, we describe a strategy, which is 
based on the Barzilai-Borwein method [2], for choosing the step size t&. Let g(-) = A*(A(-) — b) 
and g k = A* (A(X k ) — b). We perform the shrinkage iteration (|10p along the negative gradient 
direction g k of the smooth function ^||^4(X fc ) — &H2, then apply the thresholding operator T u { ) to 
accommodate the non-smooth term Hence, it is natural to choose based on the function 

\\\A(X k ) -6|| a alone. Let 

AX = X k - X k -\ Ag = g k -g k - 1 . 

The Barzilai-Borwein step provides a two-point approximation to the secant equation underlying 
quasi-Newton method, specifically, 

_ (AX,Ag) _ (AX, AX) 

Tk ~ (Ag,Ag)' °" Tk ~ (AX, Ag) 

In order to avoiding the parameter being either too small or too large, we take 

T fc = max{T min ,min{T fc ,T ma:r }}, 

where < r m j n < r max < 00 are fixed parameters. 

The idea of using the BB step to accelerate the convergence of gradient algorithms has also 
appeared in [52J. 

4.3. Algorithms 

As suggested in [T71 [3UJ, H7j , we adopt a continuation strategy to solve the regularized linear least 
squares problem ([5]). For the problem ([5]) with a target parameter fl being a moderately small 
number, we propose solving a sequence of problems ([5]) defined by an decreasing sequence [i^. 
When a new problem, associated with Hk+i, is to be solved, the approximate solution for the 
current problem with fi^ is used as the starting point. We use the parameter r\ to determine the 
rate of reduction of the consecutive fik, i-e., 

Hk+i = max^/ifej/Z), k = l,...,L—l. 

Our modified fixed point continuation iterative scheme with the Barzilai-Borwein technique for 
solving © is outlined below. 

Algorithm MFPC-BB 

Input: ► Parameters < r m j n < tq < r max < 00, /xi > p, > 0, 77 > and a tolerance e > 
Output: ► A numeric Gram matrix. 

- Set X° = 0. 

- For fi = m, . . . ,[j, L , do 
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1. Choose a step size Tfc via the BB technique such that T m j n ^ < T~max- 

2. Compute y fc = X fc -r fc ^*(^(X fc ) - 6) and a Schur decomposition of F fc = Q k A k (Q k ) T . 

3. Compute X k+1 = Q k T Tkflk (A k ) {Q k ) T '. 

- If the stop criterion is true, then return X opt . 

- end for. 

However, as shown in [3"1 \19\ [4"T]. the above algorithm may converge as 0(1/ k). Very recently, 
alternative algorithms that could speed up the performance of the gradient method FPC have been 
proposed in |19| 147] . These algorithms rely on computing the next iterate based not only on the 
previous one, but also on two or more previously computed iterates. We incorporate this new 
accelerating technique in our MFPC-BB algorithm to solve the affine constrained low-rank Gram 
matrix completion problem ([5]). The accelerated algorithm, called AFPC-BB, keeps the simplicity 
of MFPC-BB but shares the improved rate 0(1/ k 2 ) of the optimal gradient method. 

Algorithm AFPC-BB 

Input: ► Parameters < r m j n < tq < r max < oo, > p, > 0, rj > and tolerance e > 
Output: ► A numeric Gram matrix. 

- Set X° = 0. 

- For fx = m, . . . , fiL, do 

1. Choose a step size via the BB technique such that T m i n ^~max m 

2. Compute Z k = X k + ^f^(X k - X k ' x ). 

3. Compute Y k = Z k - r k A*(A(Z k ) - b) and a Schur decomposition of Y k = Q k A k {Q k ) T . 

4. Compute X k+1 = Q k T Tkf , k (A k ) {Q k ) T '. 



5. Compute t^i = 1+ ^/^ +4t k ^ 

- If the stop criterion is true, then return X op t. 

- end for. 

The following theorem shows that by performing the gradient step at the matrix Z k instead of 
at the approximate solution X k , the convergence rate of the MFPC-BB method can be accelerated 
to 0{l/k 2 ). 

Theorem 5. fM\4^ Let {X k } be the sequence generated by the AFPC-BB algorithm. Then for 
any k > 1, we have 

F(X k ) - F(X*) < " (A; + 1)2 llF , (14) 
where C is a constant, F(X) is the objective function and X* is the optimal solution of the problem 

(Sl- 



ii 



5. Numerical experiments 



In this section, we report the performance of our modified FPC algorithms for writing a real positive 
semidefinite polynomial as a sum of minimum number of squares of polynomials. In our tests, we 
generate positive semidefinite matrices W S Q nxn with rank r by randomly sampling annxr factor 
L with rational entries and setting W = LL T . After multiplying the matrix W by a monomial vector 
rrid(x) and its transpose, we obtain a positive semidefinite polynomial 



Replacing entries in W by parameters, expanding the right-hand side of the equality and matching 
coefficients of the monomials, we obtain a set of linear equations which can be written as 



where A is the linear map from S n to M. p . 

Since the SOS representation of a nonnegative polynomial is in general not unique, the so- 
lution X opt returned by MFPC-BB and AFPC-BB algorithms probably doesn't correspond to 
the constructed rational Gram matrix W. Therefore, in stead of setting relative error equal to 
||ATopt — W||.f/||W||f, which is used in [61 [301 HZ!) we choose to measure the accuracy of the 
computed solution X op t by the relative error defined by: 



The relative error also gives us a stopping criterion for the MFPC, MFPC-BB, AFPC-BB algorithms 
in our numerical experiments. We declared that the Gram matrix is approximately recovered if the 
relative error is less than a given tolerance denoted by e. 

An n x n symmetric matrix of rank r depends on d r = n(2n — r + l)/2 degrees of freedom. Let 
FR = d r /p be the ratio between the degrees of freedom in an n x n symmetric matrix of rank r and 
the number of linear constrains defined in (|15p . If FR is large (close to 1), recovering W becomes 
harder as the number of measurements is close to the degrees of freedom. Conversely, if FR is close 
to zero, recovering W becomes easier. Note that if FR > 1, there might have an infinite number of 
matrices with rank r satisfying given affine constraints. 

Throughout the experiments, we choose an initial matrix X° to be a zero matrix. For each test, 
we make an initial estimate of the value L = ||«4||| w hich is the smallest Lipschitz constant of the 
gradient of ^||*/4A" — £>|||- We set the Barzilai-Borwein parameters T max = 10/ L and T m j n = 
The thresholds 10 and 10 -3 are found after some experimentations. 

We have implemented the MFPC-BB and AFPC-BB algorithms in MATLAB, using PROPACK 
package to evaluate partial eigenvalue decompositions. All runs are conducted on a HP xw8600 
workstation with an Inter Xeon(R) 2.67GHz CPU and 3.00 GB of RAM. 

5.1. Numerical experiments on random Gram matrix completion problems 

In the first series of test, we set e = 5 x 10~ 3 and compare the performance of the MFPC, MFPC- 
BB and AFPC-BB algorithms without continuation technique to solve © for randomly generated 
matrix completion problems with moderate dimensions. In order to see the convergence behaviors 
of MFPC, MFPC-BB and AFPC-BB clearly, we compute the full Schur decompositions at each 
iteration. 

Table Q] reports the degree of freedom ratio FR, the number of iterations, and the error ()16p of 
the three algorithms MFPC, MFPC-BB, AFPC-BB. As can be seen from Table [Q on the condition 



f(x) = m d (xf ■ W • m d (x) E Q[x]. 



A(W) = b, 



(15) 



error : = 



P(^opt) - &||2 
IH|2 



(16) 
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that these three algorithms achieve similar errors, MFPC-BB provides better performance with less 
number of iterations than MFPC, while AFPC-BB outperforms the other two algorithms greatly 
in terms of the number of iterations. 



Problems 


MFPC 


MFPC-BB 


AFPC-BB 


n r p FR 


# iter error 


# iter error 


# iter error 


100 10 579 1.6494 


140 4.99e-3 


75 4.95e-3 


31 4.76e-3 


200 10 1221 1.6011 


187 4.99e-3 


105 4.97e-3 


37 4.88e-3 


500 10 5124 0.9670 


632 4.99e-3 


499 4.99e-3 


66 4.90e-3 



Table 1: Comparison of MFPC, MFPC-BB and AFPC-BB, without using continuation technique. 

In Figure 1 and Figure 2, we plot the relative error ||„4(X fc ) — 6||2/||^||2 an d approximation 
error \\X k — X pt\\p versus the iteration number of these three methods on recovering a randomly 
generated 500 x 500 Gram matrix with rank 10 respectively. We terminate these three algorithms 
when the relative error (|16p is below 5 x 10 -3 . We observe that in both cases AFPC-BB converges 
much fast than MFPC-BB and MFPC. The comparison of MFPC-BB and MFPC clearly shows 
that the Barzilai-Borwein technique is quite effective in accelerating the convergence of the MFPC 
algorithm. 



MFPC 
MFPC-BB 
- AFPC-BB 



2000 



1500 




MFPC 
MFPC-BB 
- AFPC-BB 



50 100 150 200 250 300 350 

Iteration 



50 100 150 200 250 300 350 

Iteration 



Figure 1: Relative error versus iteration Figure 2: Approximation error versus 

number for a problem with n = 500, r = iteration number for a problem with n = 

10. 500,r = 10. 

In Table [21 we report the performance of the AFPC-BB algorithm with continuation tech- 
nique on randomly generated Gram matrix completion problems. We use PRO PACK to compute 
partial eigenvalues and eigenvectors. We set the regularization parameter in problem fl3J) to be 
fi = 10 _4 ||.4*6|| and \i\ = 1/4||„4*6||. The update strategy for ^ is max(l/4/ifc_i, p) whenever the 
stopping criterion is satisfied with e = 10~ 3 . 

As indicated in the table, it takes the AFPC-BB algorithm fewer than 300 iterations on the 
average and less than 15 minutes to solve all problems in our experiments. In addition, for most 
of these problems, FR is larger than 1. Especially, FR is up to 4.5923 for the problem with 
n = 1000, r = 50. To our best knowledge, nobody has considered solving matrix completion 
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Problems 


Results 


n r p FR 


# iter time error 


100 10 579 1.6494 


76 1.48e+0 9.64e-4 


500 10 3309 1.4974 


80 2.35e+l 9.90e-4 


1000 10 10621 0.9372 


165 1.41e+2 9.95e-4 


1000 50 10621 4.5923 


120 1.10e+2 9.89e-4 


1500 10 25573 0.5848 


271 6.04e+2 9.96e-4 


1500 50 25573 2.8849 


156 4.59e+2 9.83e-4 



Table 2: Numerical results for AFPC-BB on random Gram matrix completion problems. 

problems in this situation yet. It is rather surprising that the original random Gram matrix with 
low rank can be recovered given only such a small number of affine constraints. 

5.2. Exact rational sum of squares certificates 

The numerical Gram matrix W returned by the AFPC-BB algorithm satisfies 

f(x)^m d (x) T -W-m d (x), WhO. (17) 

In order to derive an exact SOS decomposition of /, we need to start with an approximate Gram 
matrix with high accuracy. Although first-order methods are often the only practical option for 
large-scale problems, it has also been observed that the sequence {X k } computed by the AFPC-BB 
algorithm converges quite slowly to an optimal solution W*. Therefore, we apply the structure- 
preserving Gauss-Newton iterations (see [201 I21j ) to refine the Gram matrix W with low rank 
returned by the AFPC-BB algorithm: we choose a rank r which is less than or equal to the 
rank of W and compute the truncated L T DL decomposition of W to obtain an approximate SOS 
decomposition 

r 

/(^EE^") 2 - 

i=l a 

then apply standard Gauss-Newton iteration to compute Acj iCr x Q such that 

r r 

/(*) = EE c ^ x ° + Ac ^ x °) 2 + °EE Ac ^x a ) 2 )- (is) 

i=l a i=l a 

The matrix W is updated accordingly to W + AW and the iteration is stopped when the backward 
error 

9 = \\f(x)-m d {x) T ■ W-m d {x)\\ 2 (19) 

is less than the given tolerance e. If 6 remains greater than e after several Gauss-Newton iterations, 
we may increase the precision or use different r and try Gauss-Newton iterations again. After 
converting the refined matrix W into a rational matrix, we use the orthogonal projection technique 
in |20|. [21] to construct an exact rational SOS decomposition for the nonnegative polynomial /. 

It is interesting to notice that the AFPC-BB algorithm provides a low-rank Gram matrix to 
seed Gauss-Newton iterations while most of the SDP solvers GloptiPoly [18], SOSTOOLS [4"0] . 
YALMIP [29], SeDuMi 06], SDPT3 [48J and SparsePOP [50] usually return a Gram matrix with 
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maximum rank (see [22, Theorem 2.1]). For example, we consider a randomly generated Gram 
matrix completion problem with n = 200, r = 5, which is created in the same way described at the 
beginning of Section 5. The smallest 10 singular values of the numerical Gram matrix computed 
by SeDuMi are 

3.527, 2.779, 2.445, 1.369, 1.184, 0.964, 0.627, 0.101, 0.485, 0.161, 0.069. 

However, the rank of the numerical Gram matrix returned by the AFPC-BB algorithm is 14. We 
notice that by applying Gauss-Newton iterations to the low-rank Gram matrix computed by AFPC- 
BB, it is usually much easy to recover an exact SOS decomposition of the nonnegative polynomial. 

In |31j . we have used the MFPC-BB algorithm to successfully recover the exact sums of squares 
of nonnegative polynomials in |21| . 

In the following two tables, we compare the performance of the AFPC-BB algorithm and the 
SDP solver SeDuMi for recovering low rank Gram matrices from affine constraints on the same 
randomly generated examples. We also show the effectiveness of Gauss- Newton iterations run in 
Maple with Digits = 14 in refining the numerical Gram matrix. These tables report the number of 
affine constraints p, the degree of freedom ratio FR, the backward error 9, the rank of the Gram 
matrix and the running time in seconds. Table H] also shows the smallest singular value a n of 
numerical Gram matrices returned by SeDuMi. We set e = 5 x 10 -4 in the AFPC-BB algorithm, 
which is small enough to guarantee very good recoverability. 



Examples 


AFPC-BB 


Gauss-Newton iteration 


n r p FR 


rank time 


rank 9 time 


50 5 255 0.9412 


9 6.874 4.06e-l 


5 1.443e-5 4.08e+0 


100 5 579 0.8463 


9 0.860 1.75e+0 


5 1.935e-9 2.98e+l 


150 5 896 0.8259 


13 2.758 7.09e+0 


5 4.023e-8 6.28e+l 


200 5 1221 0.8108 


14 3.629 1.07e+l 


5 4.030e-5 4.69e+2 


300 5 1932 0.7712 


14 22.315 2.32e+l 


5 1.379e-9 5.61e+2 


400 5 2610 0.7624 


15 12.515 6.23e+l 


5 5.825e-5 1.22e+3 


500 5 5124 0.4859 


17 24.829 5.33e+l 


5 1.479e-5 7.92e+3 



Table 3: Exact SOS certificates via AFPC-BB and Gauss-Newton iterations. 



Examples 


SDP 


Gauss-Newton iteration 


n r p FR 


a n time 


rank 9 time 


50 5 255 0.9412 


0.701 1.03e+0 


6 3.769e-8 1.59e+l 


100 5 579 0.8463 


0.042 7.77e+0 


7 2.438e-10 6.88e+l 


150 5 896 0.8259 


0.069 1.24e+l 


7 1.883e-10 2.23e+2 


200 5 1221 0.8108 


0.069 6.58e+l 


7 4.666e-9 8.21e+2 


300 5 1932 0.7712 


0.442 2.84e+2 


7 5.679e-10 1.30e+3 


400 5 2610 0.7712 


0.114 3.94e+2 


8 9.249e-10 5.00e+3 


500 5 5124 0.4859 


0.001 2.14e+3 





Table 4: Approximate SOS certificates via SDP and Gauss-Newton iterations. 

As indicated in Table using the AFPC-BB algorithm, we can compute numerical low-rank 
Gram matrices very efficiently. Moreover, for each example, we can use Gauss-Newton iterations 
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()18p to refine the Gram matrix returned by AFPC-BB to relatively high accuracy, e.g. 10~ 5 . By 
rounding every entry of the refined matrix to the nearest integer, we can easily recover a rational 
Gram matrix with rank 5 which gives the exact SOS representation of the nonnegative polynomial. 

As indicated in Tabled for the same examples, numerical Gram matrices returned by SeDuMi 
have full rank for the given tolerance 10~ 3 , while the rank of matrices returned by AFPC-BB are 
relatively small. We seed the numerical Gram matrices returned by SeDuMi to Gauss-Newton 
iterations, the ranks of the refined matrices are always larger than 5 in order to guarantee the 
convergence of the Gauss-Newton iterations. Furthermore, we are not yet able to recover exact SOS 
decompositions even though backward errors 9 have been reduced to the order of 10 -10 . 

From Table [U it is also interesting to notice, if we decrease the degree of freedom ratio FR by 
choosing a sparse monomial vector m^(x), it is possible to recover the exact SOS representation of 
the nonnegative polynomial from the numerical low-rank Gram matrix returned by the AFPC-BB 
algorithm, without running Gauss-Newton iterations. 



Problems 


AFPC-BB 


Rational SOS 


n 


r 


P 


FR 


jf= iter 


time 


error 


time 


50 


5 


608 


0.3947 


45 


4.38e-l 


5.84e-4 


1.09e-l 


100 


5 


1167 


0.4199 


100 


1.97e+0 


8.72e-4 


3.59e-l 


150 


5 


1703 


0.4345 


217 


5.91e+0 


9.96e-4 


8.12e-l 


200 


5 


2249 


0.4402 


239 


9.11e+0 


9.99e-4 


1.50e+0 


300 


5 


3544 


0.4204 


327 


2.11e+l 


9.77e-4 


3.45e+0 


400 


10 


10078 


0.3924 


151 


2.46e+l 


9.52e-4 


1.14e+l 


500 


20 


24240 


0.4047 


142 


4.48e+l 


4.70e-4 


4.65e+l 


1000 


10 


27101 


0.3673 


436 


3.70e+2 


4.97e-4 


1.38e+2 


1000 


50 


95367 


0.5114 


395 


6.56e+2 


9.99e-5 


1.41e+3 


1500 


10 


45599 


0.3280 


554 


1.00e+3 


4.99e-4 


3.10e+2 



Table 5: Exact SOS certificates via AFPC-BB. 
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